# MEDIPS batch processing:
# extraction of QC measures to common table

library(MEDIPS)
library(BSgenome.Hsapiens.UCSC.hg19)

getwd()
# [1] "/usr/local/data/Data/NCH/MBD-Seq/MEDIP-Analysis"

# get list of files in directory
# here we're loading the stored results of the previous batch run
files <- list.files(pattern=".RData")

#####################################
# beginning of QC extraction loop
#####################################

#loop through files, process and write results to individual RData files 
for (i in files) {
	sample <- gsub("_sorted.txt.RData","",i)
	print(sample)

	# loading the RData file
	load(i)

	#####################################
	# saturation analysis
	#####################################
		
	# collating saturation analysis data for all samples
	if (which(files == i) == 1) {
		sa.ALL <- as.matrix(cbind(sa.MEDIPS.SET$numberReads, sa.MEDIPS.SET$maxEstCor[2], sa.MEDIPS.SET$maxTruCor[2]))
		colnames(sa.ALL) <- c("numberReads", "maxEstCor", "maxTruCor")
		rownames(sa.ALL)[which(files == i)] <- sample
	} else {
		sa.ALL <- rbind(sa.ALL, cbind(sa.MEDIPS.SET$numberReads, sa.MEDIPS.SET$maxEstCor[2], sa.MEDIPS.SET$maxTruCor[2]))
		rownames(sa.ALL)[which(files == i)] <- sample
	}

	#####################################
	# perform coverages analysis, 
	# i.e. coverage of CpG across genome 
	#####################################
	
	# collating coverage analysis results
	if (which(files == i) == 1) {
		ca.ALL.CpGs.covered <- as.matrix(rbind(ca.MEDIPS.SET$coveredPos[2,]))
		colnames(ca.ALL.CpGs.covered) <- c("x1","x2","x3","x4","x5","x10")
		rownames(ca.ALL.CpGs.covered)[which(files == i)] <- sample
		ca.ALL.CpGs.fraction <- as.matrix(rbind(ca.MEDIPS.SET$coveredPos[3,]))
		colnames(ca.ALL.CpGs.fraction) <- c("x1","x2","x3","x4","x5","x10")
		rownames(ca.ALL.CpGs.fraction)[which(files == i)] <- sample
	} else {
		ca.ALL.CpGs.covered <- rbind(ca.ALL.CpGs.covered, ca.MEDIPS.SET$coveredPos[2,])
		rownames(ca.ALL.CpGs.covered)[which(files == i)] <- sample
		ca.ALL.CpGs.fraction <- rbind(ca.ALL.CpGs.fraction, ca.MEDIPS.SET$coveredPos[3,])
		rownames(ca.ALL.CpGs.fraction)[which(files == i)] <- sample
	}

	#####################################
	# CpG enrichment analysis
	#####################################

	# collating enrichment analysis data for all samples
	if (which(files == i) == 1) {
		er.ALL <- as.matrix(rbind(unlist(er.MEDIPS.SET)))
		rownames(er.ALL)[which(files == i)] <- sample
	} else {
		er.ALL <- rbind(er.ALL, unlist(er.MEDIPS.SET))
		rownames(er.ALL)[which(files == i)] <- sample
	}

	rm(list = c("MEDIPS.SET", "er.MEDIPS.SET", "sa.MEDIPS.SET", "ca.MEDIPS.SET"))
	
	gc()
}

###########################################
# end of  QC extraction loop
###########################################

###########################################
# saving QC analysis results of all samples
###########################################
save(list = c("er.ALL","sa.ALL","ca.ALL.CpGs.covered","ca.ALL.CpGs.fraction"), file = "AnalysisSummaries.RData")



